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We model the expansion of an interacting atomic Bose-Einstein condensate in a disordered lat- 
tice with a nonlinear diffusion equation normally used for a variety of classical systems. We find 
approximate solutions of the diffusion equation that well reproduce the experimental observations 
for both short and asymptotic expansion times. Our study establishes a connection between the 
peculiar shape of the expanding density profiles and the microscopic nonlinear diffusion coefficients. 
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I. INTRODUCTION 

The interplay of disorder and interactions in quan- 
tum systems gives rise to a variety of interesting phe- 
nomena, ranging from glassy quantum phases to non- 
standard transport properties. In particular, interactions 
are known to be able to break the Anderson localization 
due to disorder, restoring transport in otherwise insu- 
lating systems. The prototypal systems in which such 
effect has been studied are one-dimensional disordered 
potentials, where the expansion dynamics of an initially 
localized wavepacket has been extensively investigated 
both in theory [l|-[Tl| and experiments [T3 |. There is 
now a general agreement on an anomalous diffusion pro- 
cess, where the time-exponent of the expansion is smaller 
than the one found in linear systems, e.g. < x 2 >^ t a , 
with a < 1. This subdiffusion is essentially due to the 
presence of a local diffusion coefficient that depends on 
density and therefore decreases as the wavepacket ex- 
pands. While the time-evolution of the second moment 
of the distribution observed in numerics and experiments 
agrees with microscopic models, a satisfying modeling of 
the evolution of the overall shape of the wavepacket is 
still missing. A natural question is whether this expan- 
sion process can be modeled with a nonlinear diffusion 
equation widely employed to describe related transport 
processes in classical systems [H, [13] , which contains ex- 
plicitly a density-dependent diffusion coefficient [l5l - [T8| . 
Recent numerical studies have indeed established a link 
between the nonlinear diffusion equation (NDE) and the 
asymptotic regime of subdiffusion, employing known self- 
similar solutions of the NDE [l6], [l7|- However, a com- 
parison with experimental data is not yet possible, since 
these solutions do not apply to the limited time interval 
that is possible to study in experiments, where normally 
the shape of the wavepacket changes with time [l2j . 

In this work we study this problem and derive approx- 
imate solutions of the NDE for the short-time regime ac- 
cessible in experiments. We find a relatively good agree- 
ment between the density distributions measured in the 
experiment and these solutions, and we identify a time- 
dependent exponent that links the evolution of the sec- 



ond moment to the changing shape of the distribution. 
While the present experiments lack the necessary spatial 
resolution, we find that the detailed study of such shape 
can give direct evidence of the exponent of nonlincarity 
of the local diffusion coefficient. 



II. THE DISORDERED, INTERACTING 
SYSTEM 



In the experiment we employ an ultracold cloud on 
weakly-interacting bosons in a one-dimensional optical 
lattice that mimics a truly disordered potential. More 
in detail, we realize a quasiperiodic lattice by perturbing 
a strong sinusoidal lattice with a secondary one having 
an incommensurate spacing. As described in more detail 
elsewhere [HI, , non-interacting particles in such po- 
tential can be described by the well-known Aubry- Andre 
tight-binding (single band) Hamiltonian |19| : 
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where bpbj and rij = b^bj are the standard on-site cre- 
ation, destruction and number bosonic operators, J is 
the kinetic (hopping) energy, A is the quasi-disorder 
energy and (3 is the ratio of the two lattice spacings. 
This model is known to show Anderson localization for 
A > 2 J, with an essentially energy- independent localiza- 
tion len gth £ ~ 1/ ln(A/2J), in units of the main lattice 
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spacing 

In the experiment we can realize this single-particle 
regime, and also add a controllable repulsive interaction 
between the particles, by employing potassium-39 atoms 
with a magnetically-tunable Feshbach resonance [23|, [24| . 
In presence of interaction, one needs to introduce an ad- 
ditional term in the Hamiltonian 



Him = Uy^njjrij - 1) 



(2) 



where U parameterizes the two-particles interaction en- 
ergy and Ei nt ~ Un(x,t) represents the local interaction 
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FIG. 1. (color online) Schematic representation of the 
disordered, interacting system. An initial non-interacting 
wavepacket (blue dashed line) can be decomposed into 
exponentially-localized single-particle states (red thick lines) 
of the quasiperiodic lattice (grey thin line). A weak interac- 
tion couples the states with a density-dependent rate F and 
allows an expansion of the wavepacket. Only states separated 
by one quasiperiod, l/(/3-l) , are shown for clarity. 



per particle, where n(x, t) is the time-dependent density 
distribution. Such interaction can couple distinct single- 
particle localized states, allowing for macroscopic trans- 
port. To probe the transport properties in the experi- 
ment, we initially prepare a low-temperature sample in 
the combined potential of a quasi-periodic lattice with 
A > 2 J and a harmonic trap, characterized by a Gaus- 
sian density distribution n{x). We then remove suddenly 
the trap, and let the sample expand along the lattice for 
a variable time, in presence of an additional radial con- 
finement. In Fig[T]we show a schematic representation of 
the experiment. As shown in Fig[21 we essentially observe 
no expansion if J7=0, while for finite U the distribution 
broadens and changes shape with increasing time. The 
square root of the second moment of n(x, t) increases 
with a subdiffusivc behavior of the kind 



<7(t) = vV>=^0(l+t/*0)° 1 . 



(3) 



with a characteristic exponent a in the range 0.2-0.4 for 
Eint ^ J as already discussed in ref.[l2[ (we do not con- 
sider the regime of Ei„t ^> J, where self-trapping phe- 
nomena can complicate the dynamics). While the sub- 
diffusive expansion of the width can be explained with 
heuristic models of the microscopic dynamics [HI, [2o| . 
little or no analysis is available for the evolution of the 
overall shape of n(x, t). 

Let us start discussing a simple model of the micro- 
scopic dynamics which applies to the regime of weak in- 
teractions, where one can still describe the many body- 
states as superposition of few single-particle states. A 
general expectation is that a local diffusion coefficient 
can be defined as D ~ T£ 2 , where £ is the single-particle 
localization length and T is the coupling rate of single- 
particle states by the interaction. The latter can be eval- 
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FIG. 2. (color online) a) Time evolution of the width of 
the atomic distribution for strong (blue dots) and very weak 
(black triangles) interaction. The continuous lines are fits 
with Eq.© and give a exponents of 0.28±0.02 and 0.04±0.02 
respectively; the dash-dotted line represents the asymptotic 
slope expected for normal diffusion. The initial Gaussian 
distribution (b) evolves into a flat-top distribution at longer 
times (c). The dashed lines are Gaussian fits of the tails of 
the distributions. 



uated in perturbation theory as 

27T |(z|iJ,„ t |/)| 2 

h \Ei-E f \ 



if 



(4) 



where \i) and |/) are two generic initial and final states 
(they actually represents quadruplets of single-particle 
states, because of the form of Hi nt ) and \E{ —Ef \ is their 
energy separation, which is of the order of A. Such cou- 
pling is possible only if (i\Hi nt \f) > \Ei — Et\, One can 
have two different scenarios. If Ei nt ~ A, the dominant 
couplings are within nearby states, and one finds that 
(i\Hi nt \f) is essentially Ei nt Iif, where 7j/ is an overlap 
integral of the order unity. This implies that D oc n(x,t) 2 
and, since in ID n ~ 1/<t, that D cx o~ 2 . If instead 
E int <C A, only long-distance couplings are possible, 
which tend to decrease with decreasing interaction en- 
ergy. In this case one must expect D cx n(x,t)@ cx a^ 13 , 
with (3 > 2. 

By solving the standard diffusion equation 



dn(x,t) ^d 2 n(x,t) 



dt 



dx 2 



(5) 



for the width of the distribution, da 2 (t)/dt = 2D, one 
finds a time dependence for a(t) of the form of Eq.Q, 
with a — 1/2. If the diffusion coefficient depends on the 
width itself as D cx er~' 3 , Eq.(|3]) continues to describe the 
time evolution of the width but with a time exponent 
a = 1/(2 + j3i). These simple expectations match what 
has been observed in experiments [12| and numerical cal- 
culations (e.g. Ref. [21] and references therein) for the 
evolution of a(t). 
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III. THE NONLINEAR DIFFUSION EQUATION 

Let us now turn our attention to the evolution of the 
shape of n(x,t). The idea is to start from a nonlinear 
diffusion equation (NDE) of the form 



dn(x,t) d 



Of 



dn(x, t) 
dx 



(6) 



which takes explicitly into account a density-dependent 
diffusion coefficient. The NDE is usually studied in the 
asymptotic limit, where a self-similar solution exists [26| : 



n(x, t) cx 



1/a 



for |a;| < w(t) 
for |x| > w(t) 



(7) 



The front of the diffusion w(t) has the time dependence 
w(t) oc t 1 /( 2 + a ) j and the same dependence is found for 
a{t). In our specific problem we expect an exponent a > 
2; for a = 2 it is exactly an inverted parabola. 

This self-similar solution obviously cannot reproduce 
our short-time expansion, where we clearly see a chang- 
ing shape of the distribution. As shown in Fig. 2, in the 
experiment n(x) is initially Gaussian, while later it de- 
velops a flatter top and a relatively faster decay of the 
tails. This general behavior is actually consistent with 
the picture of a density-dependent diffusion coefficient 
discussed above, which predicts a larger D at the center 
of the distribution, and a reduced one in the tails. 

We therefore look for an approximate solution of the 
NDE which can interpolate between the initial Gaussian 
distribution and the asymptotic regime. One can start 
by noting that a Gaussian can be obtained as a limit of 
a slightly different version of Eq.© 
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The conjecture is then that a solution of the form 



n(x, t) = 



AM 




b{i)x- 
w(t) 2 



l/b{t) 



\x\ < w(t)/y/b(J) 

\x\ > w{i)/JUf) 



(8) 



(9) 



might reproduce the short-time regime of the true solu- 
tion of the NDE. Here A = A(b, w) is an appropriate 
normalization coefficient, and 



6(f) = a(l - exp(-t/r)) 



(10) 



is a time-dependent exponent. To verify this conjecture, 
we solve numerically Eq. for various values of the non- 
linear diffusion exponent a, with an initial Gaussian dis- 
tribution, and we compare the calculated n(x, t) with 
the approximation above. As summarized in Figj3l we 
find that this approximation works reasonably well at all 
times, for values of the nonlinear diffusion exponent in 
the range a = 1 — 3 (see Appendix B for more details). 
In particular, the numerical n{x) is reasonably well re- 
produced by Eq.®, besides a limited deviation of the 



tails. As we will discuss later, this deviation is not an 
issue in the analysis of the experimental data, which has 
however a limited resolution. There is also a good agree- 
ment of the time evolution of the width with Eq.©, for 
an exponent a that is consistent with a = l/(a + 2). 
Finally, b(t) is reasonably well fitted by Eq.CET 
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FIG. 3. (color online) a) Time evolution of the width of the 
numerical solution of the NDE (((J for a — (black trian- 
gles), a — 2 (blue dots) and a = 3 (green squares). Lines are 
fits with Eq.© and give a exponents of 0.5, 0.24 and 0.18, 
respectively, b) Time evolution of the time-dependent expo- 
nent b, with same color and symbol scheme as above. Lines 
are fits with Eq. (|10[) . c) Fit of the numerical solution of the 
NDE for a — 2 (blue continuous) with Eq.([9} (red dashed) or 
a Gaussian (black dotted) for t=0.01; d) same as above, for 
t=0.1. 



IV. EXPERIMENTAL RESULTS 

We can now use Eq.© to fit the experimental profiles. 
FigfJ] shows the results for a set of experimental data 
with a mean initial interaction energy Ei nt = 2.3J (blue 
data in Fig. 2), with a rather good agreement. We find 
that b is close to for short expansion times (FigfJJi) and 
evolves to larger values for increasing expansion times as 
the flat-top distribution appears (FigHJa). The goodness 
of the fit for varying time can be evaluated by the coef- 
ficient of determination R 2 . Figj4}; shows the evolution 
of R 2 for both a Gaussian fit and the fit with the ap- 
proximate solution of the NDE: while the fit with the 
Gaussian gets worse as the atomic cloud expands, the fit 
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with Eq.© remains constantly good. We note that when 
the interaction energy is not strong enough to allow the 
atomic cloud expansion (black data in Fig 12]), we cannot 
appreciate a variation of the the shape of n(x), which is 
always well fitted by a Gaussian. 
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FIG. 4. (color online) a,b) Fit of the experimental profiles 
(blue continuous line) with Eq.© (red dashed line). The 
initial Gaussian distribution (a) is well fitted with b « 0, while 
at longer times (b) the exponent increases, c) Coefficient of 
determination R 2 for a fit with Eq.([9)l (red triangles) and for 
a Gaussian fit (black dots). The lines are a guide to the eye. 

We can now compare the evolution of the exponents 
b for the experiment and the numerical solution of the 
NDE. A direct comparison can be obtained by studying 
the evolution of 6 as a function of the rescaled width 
a(t)/ao, as shown for example in FigJS] This allows to 
get rid of the different time and width scales in the ex- 
periment and in the simulations. The experimental data 
refer to the same set of FigHl for which we measured 
a = 0.28 ±0.02, it is therefore natural to compare it with 
the solution of the NDE for a = 2. One can note a quali- 
tatively similar behavior of theory and experiment, with 
an initial b « that increases towards an asymptotic 
value. However, the asymptotic value for the experiment 
is not the expected one, and the overall evolution is ap- 
parently slower. 

We attribute this discrepancy to the finite spatial reso- 
lution with which the experimental n(x) is detected. Ac- 
tually, for our imaging system we have a Gaussian point 
spread function with a width of <7r=10/im, which is there- 
fore comparable to the initial width cto- The expected 
effect of such finite resolution is indeed a smoothing of 
the steep decay of the tails of n(x), and therefore a de- 
crease of the measured exponent b. A comparison of the 
numerical and experimental profiles in Figs. 3-4 confirms 
this argument. 

A more quantitative comparison can be made by 
properly taking the finite resolution into account. To 
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FIG. 5. (color online) Evolution of the exponent b with the 
rescaled width obtained by fitting n(x, t) with the approxi- 
mate solution of the NDE: experiment (blue triangles) numer- 
ical solution of the NDE with a — 2 (black line), numerical 
solution of the DNLSE (red circles). 



strengthen our analysis, we also performed a numer- 
ical simulation of the expansion by employing a one- 
dimensional discrete non-linear Schrodinger equation 
(DNLSE) that is known to reproduce the evolution of 
our type of disordered system in the regime of weak in- 
teraction. One example of the numerical n(t) for a long 
expansion time is shown in Fig[6] In absence of a finite 
spatial resolution (FigJB^), one notes rather steep tails 
that can indeed be fitted with the approximate solution 
of the NDE with an exponent b w 2. Actually, the full 
evolution of b(t) for the solution of the DNLSE, which 
is also reported in FigjSJ shows a rather good agreement 
with the solution of the NDE at all times. When instead 
the distribution is convolved with the calculated Gaus- 
sian transfer function (Figj6j)), one can observe a clear 
smoothing of the tails, leading to a substantial reduction 
of b. 



■ a) I b= 1.7 ±0.2 
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FIG. 6. (color online) Effect of the finite imaging resolution 
on the numerical simulations with the DNLSE. a) n(x) after 
a large expansion time (black continuous line), in absence of 
a finite resolution, fitted by Eq.@ (red dashed line), b) the 
same n(x) after convolution with Gaussian transfer function 
(gray dotted line). 

A direct comparison of the experiment with the NDE 
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can therefore be made only by properly taking into ac- 
count such finite resolution also in the numerical solution 
of the NDE. Fig 17] compares the exponent b from the ex- 
periment with those fitted from the numerical solutions 
of both NDE and DNLSE, convolved with the Gaussian 
point-spread function. The evolution of the NDE expo- 
nent is now slower, and rather close to the experimental 
one. Clearly, the asymptotic value b = a can be reached 
only if the width of the distribution becomes much larger 
than a/. 



is characterized by steep tails. We can see a qualitative 
agreement on the evolution of the shape of the distri- 
bution, which confirms the hypothesis of a microscopic 
density-dependent diffusion coefficient. The finite exper- 
imental spatial resolution did not so far allow us to ver- 
ify the expected relation between the spatial and time 
exponents. This might however be explored in future 
experiments having state-of-the-art spatial resolution. 




FIG. 7. (color online) Evolution of the exponent b with the 
rescaled width obtained by fitting n(x, t), taking into account 
the finite spatial resolution: experiment (blue triangles), nu- 
merical solution of the NDE with a — 2 (black line), numerical 
solution of the DNLSE (red circles). 

The specific set of experimental data we have discussed 
so far is just one example. We have actually found sim- 
ilar results for other values of the initial interaction en- 
ergy in the range Ei nt =0.5-3J. For weaker interactions, 
although a small expansion can be detected in the ex- 
periment (see for example the data in Fig. 2), the shape 
stays essentially Gaussian up to the longest observation 
time. We speculate that other effects existing in the ex- 
perimental setups, such as a weak time-dependent noise, 
might be responsible for this observation [271] . 

V. CONCLUSIONS 



VI. APPENDICES 

A. Experimental methods and parameters 

The quasiperiodic potential is created by perturbing a 
deep optical lattice with a weaker lattice of incommen- 
surate wavelength: V(x)—V\ cos 2 (fcix) + V2Cos 2 (A)2a; + 
4>). Here ki=2ir/Ai are the wavevectors of the lattices, 
with Ai = 1064.4 nm and A2=859.6 nm, giving a ratio 
(3— 1.238..., 4> is the relative phase between the two lat- 
tices. The main lattice fixes the lattice spacing, d=Ai/2, 
and the tunneling energy J. The quasi-disorder strength, 
A, scales linearly with the secondary lattice strength Va 

The atomic sample consists in a Bose-Einstein con- 
densate of 39 K atoms in their ground state, whose s- 
wave scattering length a s can be adjusted from about 
zero to large positive values thanks to a broad magnetic 
Feshbach resonance [23|, [24| . The condensate is initially 
produced in a crossed optical trap at a s = 280a , and 
contains about 4 x 10 4 atoms. A quasiperiodic lattice 
with A « 3J is then slowly added. The radial confine- 
ment induced by the optical trap and the lattice beams is 
u) r s=! 2ir x 80 Hz while the axial one is cu ax « 2ir x 70 Hz. 

At a given time, t = 0, the optical trap is suddenly 
switched off and the atoms are let free to expand along 
the lattice, in presence of a radial confinement of u r ~ 
2ir x 50 Hz given by the radial profile of the lattice beam. 
At the same time, A and a are tuned to their final values 
within 10 ms, and kept there for the rest of the evolution. 
The interaction parameter is set by the scattering length 
as 



In conclusion, we have compared the evolution of the 
density distribution during the subdiffusive expansion of 
interacting atoms in disorder with the solution of a gen- 
eral non-linear diffusion equation of the form of Eq. (j6j) . 
In this equation the local diffusion coefficient is consid- 
ered to be proportional to the density to some power a. 
To account for the relatively small timescales available 
in experiments, we have built an approximate solution 
of the nonlinear diffusion equation that provides a suffi- 
ciently accurate interpolation between the initial Gaus- 
sian wavepacket and the asymptotic distribution, which 



U = / (p 4 d 3 x, (11) 

m J 

where (p is the 3D single-particle Wannier wavefunction. 
A maximum Ei nt ~ J can be realized, since an increas- 
ing repulsion tends to broaden the system radially, thus 
reducing its density. 

The atomic density distribution is detected via absorp- 
tion imaging, and integrated along the radial direction to 
obtain the one-dimensional profiles n(x). 
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B. Approximate solution of the nonlinear diffusion 
equation 

The complete expression for the approximation to the 
solutions of the NDE, normalized to unity, is 



n(x, t) 



6 3 / 2 r(l/6 + 3/2) / k 2 \iA 



y/^wT(l/b) 



1 - 



(12) 



for | a: | < w/Vb, and zero otherwise. This expression 
provides an overall better fit of the numerical solutions 
of the NDE in all time regimes than either a Gaussian 
or the asymptotic solution of the NDE. Fig. 8 shows for 
example the coefficient of determination for the specific 
case of a = 2. 



Pi 1.04 




FIG. 8. (color online) Coefficient of determination R 2 of the 
fit of the numerical solutions of the NDE with Eq. (|12|l (blue 
triangles), with a Gaussian (black squares) or with Eq.© (red 
dots). 



C. Numerical simulations 

The DNLSE studied in the numerical simulations re- 
produces the experimental system in the limit of negligi- 
ble population of the radial degrees of freedom: 



-(^+i(t)+v,-i(t))+ 



+A/Jsin(27r/3j + 4>)\^{t)\ 2 + 7r%lVj ■ (13) 

Here ipj(t) are the coefficients of the wave function in 
the Wannier basis, normalized in such a way that their 
squared modulus corresponds to the atom density on 
the j-th site of the lattice. The mean field interaction 
strength is given by 7. The relation between Ei nt and 7 
is approximately Ei nt « 2Jj/n s , where n s is the mean 
number of sites occupied by the atomic distribution. The 
initial condition of the simulations is a Gaussian distribu- 
tion as in the experiment. For each expansion time, n(x) 
is obtained by averaging the profiles resulting from 100 
different realizations of the quasiperiodic potential with 
randomly varied phase <j> in the range [0, 2ir). The spe- 
cific results of Figs. 5-7 where obtained with parameters 
7 = 40 and A/J = 2.5 
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